/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of LINAL (LINear ALgebra)
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "linal/common.hxx"
#include "linal/util.hxx"

#include "linal/format.hxx"

#include "linal/device/host/operation.hxx"

namespace linal {
  namespace host {
    namespace internal {

      int apply_diag_scal       ( int side, int conj, int inv, int uplo, FLA_Obj &A ) { 

        FLA_Obj ATL,   ATR,      A00,  a01,      A02,
                ABL,   ABR,      a10,  alpha11,  a12,
                                 A20,  a21,      A22;

        FLA_Obj beta;

        FLA_Obj_create( FLA_Obj_datatype(A), 1, 1, 0, 0, &beta );

        FLA_Part_2x2( A, &ATL, &ATR,
                         &ABL, &ABR, 0, 0, FLA_TL );

        while ( FLA_Obj_length( ATL ) < FLA_Obj_length( A ) ){

          FLA_Repart_2x2_to_3x3( ATL, ATR,     &A00, &a01,     &A02,
                                               &a10, &alpha11, &a12,
                                 ABL, ABR,     &A20, &a21,     &A22,
                                 1, 1, FLA_BR );
          //-----------------------------------------------------------

          FLA_Set( alpha11, beta );

          switch (inv) {
          case FLA_INVERSE:
            FLA_Invert( conj, beta );
            break;
          }

          switch (side) {
          case FLA_LEFT:
            switch (uplo) {
            case FLA_LOWER_TRIANGULAR:
              FLA_Scal( beta, a10 );
              break;
            case FLA_UPPER_TRIANGULAR: 
              FLA_Scal( beta, a12 );
              break;
            }
            break;
          case FLA_RIGHT:
            switch (uplo){
            case FLA_LOWER_TRIANGULAR: 
              FLA_Scal( beta, a21 );
              break;
            case FLA_UPPER_TRIANGULAR: 
              FLA_Scal( beta, a01 );
              break;
            }
            break;
          }
          
          //-----------------------------------------------------------
          FLA_Cont_with_3x3_to_2x2( &ATL,  &ATR,  A00, a01,     A02,
                                                  a10, alpha11, a12,
                                    &ABL, &ABR,   A20, a21,     A22,
                                    FLA_TL );
        }

        FLA_Obj_free( &beta );

        return true; 
      }

      int apply_diag_scal       ( int side, int conj, int uplo, FLA_Obj &A ) { 
        return apply_diag_scal( side, conj, FLA_NO_INVERSE, uplo, A );
      }
      int apply_diag_inv_scal   ( int side, int conj, int uplo, FLA_Obj &A ) { 
        return apply_diag_scal( side, conj, FLA_INVERSE, uplo, A );
      }

      int apply_diag_scal       ( int side, int conj, int inv, FLA_Obj &A, FLA_Obj &B ) {

        Flat_ AA(A), BB(B), b, alpha, gamma;
        int is_mat = (AA.get_n() > 1), m = BB.get_m(), n = BB.get_n();

        gamma.create( AA.get_datatype(), 1, 1 );

#pragma unroll
        for (int i=0;i<AA.get_m();++i) {
          
          AA.extract( alpha, 1, 1, i, (is_mat ? i:0) );
          FLA_Set( ~alpha, ~gamma );

          switch (inv) {
          case FLA_INVERSE: FLA_Invert( conj, ~gamma );             break;
          }

          switch (side) {
          case FLA_LEFT:    BB.extract( b, 1, n, i, 0 );            break;
          case FLA_RIGHT:   BB.extract( b, m, 1, 0, i );            break;
          }
          FLA_Scal( ~gamma, ~b );
        }

        gamma.free();

        return true; 
      }
      int apply_diag_scal       ( int side, int conj, FLA_Obj &A, FLA_Obj &B ) {
        return apply_diag_scal( side, conj, FLA_NO_INVERSE, A, B ); 
      }
      int apply_diag_inv_scal   ( int side, int conj, FLA_Obj &A, FLA_Obj &B ) {
        return apply_diag_scal( side, conj, FLA_INVERSE, A, B ); 
      }
    }
  } 
}
